Install any packages you might be missing:
library(spatstat)
library(splancs)
# tree data
data(bramblecanes)
data(finpines)
data(spruces)
#get info
print(finpines)
## Marked planar point pattern: 126 points
## Mark variables: diameter, height
## window: rectangle = [-5, 5] x [-8, 2] metres
intensity(finpines)
## [1] 1.26
summary(finpines)
## Marked planar point pattern: 126 points
## Average intensity 1.26 points per square metre
##
## Coordinates are given to 13 decimal places
##
## Mark variables: diameter, height
## Summary:
## diameter height
## Min. :0.000 Min. :0.800
## 1st Qu.:1.000 1st Qu.:1.825
## Median :2.000 Median :2.850
## Mean :2.532 Mean :2.828
## 3rd Qu.:3.000 3rd Qu.:3.600
## Max. :7.000 Max. :5.400
##
## Window: rectangle = [-5, 5] x [-8, 2] metres
## Window area = 100 square metres
## Unit of length: 1 metre
#creates point pattern data format from the x and y coordinates
finpines.pts<-as.points(finpines$x,finpines$y)
bramblecanes.pts <-as.points(bramblecanes$x,bramblecanes$y)
spruces.pts <-as.points(spruces$x,spruces$y)
dev.off()
## null device
## 1
#divide into quadrants and count number of points in the quadrant
Q <- quadratcount(finpines, nx = 4, ny = 3)
Q
## x
## y [-5,-2.5) [-2.5,0) [0,2.5) [2.5,5]
## [-1.33,2] 10 8 11 21
## [-4.67,-1.33) 8 7 9 3
## [-8,-4.67) 12 11 22 4
plot(finpines.pts,pch=19, main="Finpines", xlab="Easting",ylab="Northing")
plot(Q, add = TRUE, cex = 2)
# Finpines data
sp_point <- matrix(NA, nrow=length(finpines$x),ncol=2)
sp_point[,1] <- finpines$x
sp_point[,2] <- finpines$y
colnames(sp_point) <- c("x","y")
plot(x=sp_point[,1],y=sp_point[,2],main="Finpines Data", xlab="Easting",ylab="Northing",cex=.5)
intensity(finpines)
## [1] 1.26
summary(finpines)
## Marked planar point pattern: 126 points
## Average intensity 1.26 points per square metre
##
## Coordinates are given to 13 decimal places
##
## Mark variables: diameter, height
## Summary:
## diameter height
## Min. :0.000 Min. :0.800
## 1st Qu.:1.000 1st Qu.:1.825
## Median :2.000 Median :2.850
## Mean :2.532 Mean :2.828
## 3rd Qu.:3.000 3rd Qu.:3.600
## Max. :7.000 Max. :5.400
##
## Window: rectangle = [-5, 5] x [-8, 2] metres
## Window area = 100 square metres
## Unit of length: 1 metre
# Random points
u.x <- runif(n=nrow(sp_point), min=bbox(sp_point)[1,1], max=bbox(sp_point)[1,2])
u.y <- runif(n=nrow(sp_point), min=bbox(sp_point)[2,1], max=bbox(sp_point)[2,2])
# Regular points
r.x <- seq(from=min(sp_point[,1]),to=max(sp_point[,1]),length=sqrt(nrow(sp_point)))
r.y <- seq(from=min(sp_point[,2]),to=max(sp_point[,2]),length=sqrt(nrow(sp_point)))
r.x <- jitter(rep(r.x,length(r.x)),.001)
r.y <- jitter(rep(r.y,each=length(r.y)),.001)
# Plot Finpines, uniform, and regular points
par(mfrow=c(1,3),mar=c(4,4,1.5,0.5))
plot(x=sp_point[,1],y=sp_point[,2],main="Finpines Data", xlab="Easting",ylab="Northing",cex=.5,pch=19)
plot(x=u.x,y=u.y,main="Random Points", xlab="Easting",ylab="Northing",cex=.5,pch=19)
plot(x=r.x,y=r.y,main="Regular Points", xlab="Easting",ylab="Northing",cex=.5,pch=19)
# CSR from random uniform point process
## Generates a CSR point process conditionally on a fixed number N of points
## Let N = 100
par(mfrow=c(2,2),mar=c(0.4,0.4,1,0.4))
pp.100.unif<-runifpoint(100,win=owin(c(0,1),c(0,1)))
pp.100.unif
## Planar point pattern: 100 points
## window: rectangle = [0, 1] x [0, 1] units
intensity(pp.100.unif)
## [1] 100
plot(pp.100.unif, main="")
title('100 uniform points')
# repeat a few times to see different realizations of the process
# HPP with intensity 100 in the unit square
pp.100.unit <- rpoispp(100,win=owin(c(0,1),c(0,1)))
pp.100.unit
## Planar point pattern: 107 points
## window: rectangle = [0, 1] x [0, 1] units
intensity(pp.100.unit)
## [1] 107
plot(pp.100.unit,main="")
title('HPP intensity = 100, unit square')
# repeat a few times to see different realizations of the process (including different numbers of points)
# HPP with intensity 1 in a 10 x 10 square
pp.1.10 <- rpoispp(1, win=owin(c(0,10),c(0,10)))
par(mfrow=c(1,2))
plot(pp.100.unit,main="")
title('intensity = 100, unit square')
plot(pp.1.10,main="")
title('intensity = 1, 10 x 10 square')
# HPP with intensity 200 in the unit square
par(mfrow=c(1,1))
pp.200.unit <- rpoispp(200)
plot(pp.200.unit)
# HPP with intensity 100 in the unit square
side.square <-1
pp.100.unit <- rpoispp(10000,win=owin(c(0,side.square),c(0,side.square)))
points.tot <- npoints(pp.100.unit)
plot(pp.100.unit,main="")
# inscribe a circle within the square
radius.circle <- side.square/2
win.circle <- disc(radius.circle, c(radius.circle,radius.circle))
pp.incircle <- pp.100.unit[win.circle]
points.incircle <- npoints(pp.incircle)
summary(pp.incircle)
## Planar point pattern: 7942 points
## Average intensity 10116.13 points per square unit
##
## Coordinates are given to 8 decimal places
##
## Window: polygonal boundary
## single connected closed polygon with 128 vertices
## enclosing rectangle: [0, 1] x [0, 1] units
## Window area = 0.785083 square units
## Fraction of frame area: 0.785
plot(pp.incircle)
# calculate pi using ratio of areas and number of points
piEst <- ((side.square^2)*points.incircle) / ((radius.circle^2)*points.tot)
piEst
## [1] 3.13913
Function to simulate pi over many iterations
pi.simulation2 = function(iter, lambda, side.square){
piEst.list = rep(0,iter)
points.incircle.list = rep(0,iter)
# inscribe a circle within the square
radius.circle <- side.square/2
win.circle <- disc(radius.circle, c(radius.circle,radius.circle))
for(i in 1:iter){
# HPP with intensity 100 in the unit square
pp.square <- rpoispp(lambda,win=owin(c(0,side.square),c(0,side.square)))
points.tot <- npoints(pp.square)
# find points within the circle
pp.incircle <- pp.square[win.circle]
points.incircle <- npoints(pp.incircle)
# calculate pi using ratio of areas and number of points
pihat <- ((side.square^2)*points.incircle) / ((radius.circle^2)*points.tot)
piEst.list[i] = pihat
}
return( piEst.list)
}
iter = 1000
lambda = 1000
side.square = .5
est = pi.simulation2(iter,lambda, side.square)
mean(est)
## [1] 3.144389
var(est)
## [1] 0.01097336
# we plot the histogram and see it appears to be binomial. why?
hist(est,nclass = 15)
# create polygon (square) area for finpines data
minfx <- min(finpines$x)
maxfx <- max(finpines$x)
minfy <- min(finpines$y)
maxfy <- max(finpines$y)
polyfin <- as.points(c(minfx,maxfx,maxfx,minfx),c(minfy,minfy,maxfy,maxfy))
UL.khat.fin <- Kenv.csr(length(finpines$x), polyfin, nsim=99, seq(0.05,3.5,0.15))
## Doing simulation 1
## Doing simulation 2
## Doing simulation 3
## Doing simulation 4
## Doing simulation 5
## Doing simulation 6
## Doing simulation 7
## Doing simulation 8
## Doing simulation 9
## Doing simulation 10
## Doing simulation 11
## Doing simulation 12
## Doing simulation 13
## Doing simulation 14
## Doing simulation 15
## Doing simulation 16
## Doing simulation 17
## Doing simulation 18
## Doing simulation 19
## Doing simulation 20
## Doing simulation 21
## Doing simulation 22
## Doing simulation 23
## Doing simulation 24
## Doing simulation 25
## Doing simulation 26
## Doing simulation 27
## Doing simulation 28
## Doing simulation 29
## Doing simulation 30
## Doing simulation 31
## Doing simulation 32
## Doing simulation 33
## Doing simulation 34
## Doing simulation 35
## Doing simulation 36
## Doing simulation 37
## Doing simulation 38
## Doing simulation 39
## Doing simulation 40
## Doing simulation 41
## Doing simulation 42
## Doing simulation 43
## Doing simulation 44
## Doing simulation 45
## Doing simulation 46
## Doing simulation 47
## Doing simulation 48
## Doing simulation 49
## Doing simulation 50
## Doing simulation 51
## Doing simulation 52
## Doing simulation 53
## Doing simulation 54
## Doing simulation 55
## Doing simulation 56
## Doing simulation 57
## Doing simulation 58
## Doing simulation 59
## Doing simulation 60
## Doing simulation 61
## Doing simulation 62
## Doing simulation 63
## Doing simulation 64
## Doing simulation 65
## Doing simulation 66
## Doing simulation 67
## Doing simulation 68
## Doing simulation 69
## Doing simulation 70
## Doing simulation 71
## Doing simulation 72
## Doing simulation 73
## Doing simulation 74
## Doing simulation 75
## Doing simulation 76
## Doing simulation 77
## Doing simulation 78
## Doing simulation 79
## Doing simulation 80
## Doing simulation 81
## Doing simulation 82
## Doing simulation 83
## Doing simulation 84
## Doing simulation 85
## Doing simulation 86
## Doing simulation 87
## Doing simulation 88
## Doing simulation 89
## Doing simulation 90
## Doing simulation 91
## Doing simulation 92
## Doing simulation 93
## Doing simulation 94
## Doing simulation 95
## Doing simulation 96
## Doing simulation 97
## Doing simulation 98
## Doing simulation 99
khat.fin<-khat(finpines.pts, polyfin, seq(0.05,3.5,0.15))
# plot of Khat-pi t^2 from data
plot(seq(0.05,3.5,0.15), khat.fin-pi*seq(0.05,3.5,0.15)^2, type="l", xlab="Distance", ylab="Estimated K-pi h^2")
# plot upper bound
lines(seq(0.05,3.5,0.15), UL.khat.fin$upper-pi*seq(0.05,3.5,0.15)^2, lty=2)
# plot lower bound
lines(seq(0.05,3.5,0.15), UL.khat.fin$lower-pi*seq(0.05,3.5,0.15)^2, lty=2)
# L functions
l<-function(k,h)
{sqrt(k/pi)-h}
# plot Lhat from data
plot(seq(0.05,3.5,0.15), l(khat.fin, seq(0.05,3.5,0.15)), type="l", xlab="Distance", ylab="Estimated L",ylim=c(-.1,.2))
# plot upper bound of Lhat
lines(seq(0.05,3.5,0.15), l(UL.khat.fin$upper,seq(0.05,3.5,0.15)), lty=2)
# plot lower bound of Lhat
lines(seq(0.05,3.5,0.15), l(UL.khat.fin$lower,seq(0.05,3.5,0.15)), lty=2)
K<-Kest(finpines)
plot(K)
plot(K, cbind(r, sqrt(iso/pi)) ~ r)
plot(K, cbind(trans,iso,border) - theo ~ r)
plot(envelope(finpines,Kest))
## Generating 99 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
## 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
## 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99.
##
## Done.
# G test: Nearest neighbor distribution test
Gtest<-Gest(finpines)
plot(envelope(finpines,Gest))
## Generating 99 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
## 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
## 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99.
##
## Done.
# F test: Point to nearest event distribution test
Ftest<-Fest(finpines)
plot(envelope(finpines,Fest))
## Generating 99 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
## 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
## 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99.
##
## Done.
# J test: J(r) = (1-G(r)/(1-F(r)))
Jtest<-Jest(finpines)
plot(envelope(finpines,Jest))
## Generating 99 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
## 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
## 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99.
##
## Done.
lamb1 <- function(x,y) {100 * exp(10*x-5*y)}
Ipp1 <- rpoispp(lamb1, 100)
Ipp1.p<- as.points(Ipp1$x,Ipp1$y)
plot(Ipp1.p,main="IPP, intensity(x,y)=100*exp(10x-5y)",pch=20,xlab='x',ylab='y')
lamb2 <- function(x,y) {100 * exp(-10*x+5*y)}
Ipp2 <- rpoispp(lamb2, 100)
Ipp2.p<- as.points(Ipp2$x,Ipp2$y)
plot(Ipp2.p,main="IPP, intensity(x,y)=100*exp(-10x+5y)",pch=20,xlab='x',ylab='y');
# can save a step and write function within rpoispp
Ipp3 <- rpoispp(function(x,y) {100 * exp(-3*x)}, 100)
plot(Ipp3,pch=20,xlab='x',ylab='y')
#Generate a square polygon
spoly<-as.points(c(0,1,1,0),c(0,0,1,1)) #polygon
#intensity of parents is 15, average number of offspring is 7, and variance of distance traveled by child (away from parent) in 0.0025
pcp1<-pcp.sim(15,7,0.00025,spoly)
plot(pcp1,main="PCP, (P,O,Spread)=(25,7,.00025)",pch=20,xlab='x',ylab='y')
#intensity of parents is 25, average number of offspring is 4, and variance of distance traveled by child in 0.005
pcp2 <-pcp.sim(25,4,0.005,spoly)
plot(pcp2,main="PCP, (P,O,Spread)=(25,4,.005)",pch=20,xlab='x',ylab='y')
# Using K-function for Clustered Process
# Simulate homogeneous pcp (same as above)
pcp2 <-pcp.sim(25,4,0.005,spoly)
# need to define polygonal area
spoly<-as.points(c(0,1,1,0),c(0,0,1,1))
# estimate
pcp.fit<-pcp(pcp2,spoly,h0=.3,n.int=40)
#simulate from the theoretical process with parameter estimates from above
plot(pcp2)
m.h <- npts(pcp2)/(areapl(spoly)*pcp.fit$par[2])
sims.h <- pcp.sim(pcp.fit$par[2], m.h, pcp.fit$par[1], spoly)
pointmap(as.points(sims.h), add=TRUE, col="red")
# L-function
l<-function(k,s)
{sqrt(k/pi)-s}
# K function for PCP
r <- seq(0.01,1,by=.2)
K.env <- Kenv.pcp(pcp.fit$par[2], m.h, pcp.fit$par[1], spoly, nsim=20, r=r)
L.env<- lapply(K.env,FUN=l,r)
limits <- range(unlist(L.env))
plot(r, sqrt(khat(as.points(pcp2),spoly,r)/pi)-r, ylim=limits, main="L function with simulation envelopes and average", type="l", xlab="distance", ylab="")
lines(r, L.env$lower, lty=5)
lines(r, L.env$upper, lty=5)
lines(r, L.env$ave, lty=6)
abline(h=0)
data(cardiff)
cardiff.pts <-as.points(cardiff$x,cardiff$y)
polymap(cardiff$poly)
pointmap(as.points(cardiff), add=TRUE)
title("Locations of homes of 168 juvenile offenders")
pcp.fit.cardiff <- pcp(as.points(cardiff), cardiff$poly, h0=20, n.int=30)
pcp.fit.cardiff
## $par
## s2 rho
## 6.32912666 0.01153204
##
## $value
## [1] 0.02499235
##
## $counts
## function gradient
## 77 NA
##
## $convergence
## [1] 0
##
## $message
## NULL
m <- npts(as.points(cardiff))/(areapl(cardiff$poly)*pcp.fit.cardiff$par[2])
r <- seq(2,30,by=2)
K.env <- Kenv.pcp(pcp.fit.cardiff$par[2], m, pcp.fit.cardiff$par[1], cardiff$poly,
nsim=20, r=r)
L.env<- lapply(K.env,FUN=l,r)
limits <- range(unlist(L.env))
plot(r, sqrt(khat(as.points(cardiff),cardiff$poly,r)/pi)-r, ylim=limits,
main="L function with simulation envelopes and average", type="l",
xlab="distance", ylab="")
lines(r, L.env$lower, lty=5)
lines(r, L.env$upper, lty=5)
lines(r, L.env$ave, lty=6)
abline(h=0)